
;;=================================================================
;; Load the GRS thang...
angsep = 10.
;; Get galactic params
defsysv, '!MW', exists = exists
IF ~ exists THEN galactic_params 
d = (n_elements(dvec) NE 0) ? dvec : $
    dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART

restore,'./ancillary/grs_kda_data.sav',/ver
restore,'./emaf_paper/bgps_emaf.sav',/ver
n = n_elements(emaf)

restore,'./ancillary/grs_kda_data.sav',/ver
restore,'./emaf_paper/bgps_emaf.sav',/ver
n = n_elements(emaf)

;; Create arrays used for GRS identification
dsep = fltarr(n)
vsep = fltarr(n)
grsi = lonarr(n)

;; Loop through EMAF to check for GRS co-incidence
message,'Checking for GRS co-incidence...',/inf
FOR i=0L, n-1 DO BEGIN
   gcirc, 2, emaf[i].l, emaf[i].b, grs.glon, grs.glat, dis
   mind = min(dis,gind)
   vsep[i] = abs(emaf[i].vlsr - grs[gind].vlsr)
   
   ;; A DSEP =-1 means there is no coincidence for this object
   dsep[i] = (vsep[i] LE grs[gind].delv*1.0) ? mind/60. : -1
   IF dsep[i] GT angsep THEN dsep[i] = -1
   grsi[i] = gind
ENDFOR

emaf.grs = dsep GE 0.
message,'There are '+string(total(emaf.grs),format="(I0)")+$
        ' GRS-associated sources.',/inf

;; Load emaf
emaf.grs_clemdist = grs[grsi].dis
emaf.grs_v        = grs[grsi].vlsr
emaf.grs_l        = grs[grsi].glon
emaf.grs_b        = grs[grsi].glat

;; SELECT ONLY GRS-ASSOCIATED SOURCES
ind = where(emaf.grs, n)
emaf = emaf[ind]

;; Derive GRS distance based on Reid (2009) rotation curve
;; Read in GRS distance correction file
readcol,'./emaf_paper/corrected_grs_distance.txt',bad_gd_cnum,new_gd,$
        format='I,F'
bgdind = WHERE_ARRAY(bad_gd_cnum,emaf.cnum,nbgd)
ngdind = WHERE_ARRAY(emaf.cnum,bad_gd_cnum,nngd)
emaf[bgdind].grs_clemdist = new_gd[ngdind]
grs_dtan = !MW.R0 * cos(emaf.grs_l * !dtor) / cos(emaf.grs_b * !dtor)/1.d3
grs_far = emaf.grs_clemdist GE grs_dtan
print,'N GRS at far KDIST: ',long(total(grs_far))


;; Calculate Reid distances for GRS sources
FOR i=0L,n-1 DO BEGIN
   IF (i+1) MOD 20 EQ 0 THEN $
      print,'Processing object '+string(i+1,format="(I3)")+$
            ' of '+string(n,format="(I0)")
   
   ;; Calculate grs_dist based on GRS KDA resolution
   emaf[i].grs_dist=(grs_far[i]) ? $
                    KDIST(emaf[i].grs_l,emaf[i].grs_b,emaf[i].grs_v,/FAR)/1.d3:$
                    KDIST(emaf[i].grs_l,emaf[i].grs_b,emaf[i].grs_v,/NEAR)/1.d3
ENDFOR

;;=================================================================
;; Back to your regularly-scheduled computations...

nup      = 250
upsilon  = dindgen(nup)*.01d + 0.01d
grs_dml  = dblarr(nup)
grs_dbar = dblarr(nup)
dml_c1   = dblarr(nup)
dml_c15  = dblarr(nup)
dml_c2   = dblarr(nup)
dbar_c1  = dblarr(nup)
dbar_c15 = dblarr(nup)
dbar_c2  = dblarr(nup)


indc1  = WHERE( emaf.c GE 0.10,  nc1 )
indc15 = WHERE( emaf.c GE 0.15, nc15 )
indc2  = WHERE( emaf.c GE 0.20,  nc2 )

geq = cgSymbol('geq',/ps)

FOR i=0,nup-1 DO BEGIN
   
   message,'Processing for upsilon = '+string(upsilon[i],format="(F0.2)"),/inf
   
   fn = './emaf_paper/save_files/morph_distances_'+$
        string(upsilon[i],format="(F0.2)")+'.sav'
   IF ~ FILE_TEST(fn) THEN $
      run_morph_matching,/drivemorph,upsilon=upsilon[i],/nops,/silent
   ;; fn = './emaf_paper/morph_distances.sav'
   restore,fn,/ver
   
   dml = dml[ind,*]
   dbar = dbar[ind,*]
   
   dmlind  = where( abs(dml[*,0]  - emaf.grs_dist) LE 1., nml)
   dbarind = where( abs(dbar[*,0] - emaf.grs_dist) LE 1., nbar)
   grs_dml[i]  = double(nml)  / double(n)
   grs_dbar[i] = double(nbar) / double(n)
   
   
   dmlind  = where( abs(dml[indc1,0]  - emaf[indc1].grs_dist) LE 1., nml)
   dbarind = where( abs(dbar[indc1,0] - emaf[indc1].grs_dist) LE 1., nbar)
   dml_c1[i]  = double(nml)  / double(nc1)
   dbar_c1[i] = double(nbar) / double(nc1)
   
   dmlind  = where( abs(dml[indc15,0]  - emaf[indc15].grs_dist) LE 1., nml)
   dbarind = where( abs(dbar[indc15,0] - emaf[indc15].grs_dist) LE 1., nbar)
   dml_c15[i]  = double(nml)  / double(nc15)
   dbar_c15[i] = double(nbar) / double(nc15)
   
   dmlind  = where( abs(dml[indc2,0]  - emaf[indc2].grs_dist) LE 1., nml)
   dbarind = where( abs(dbar[indc2,0] - emaf[indc2].grs_dist) LE 1., nbar)
   dml_c2[i]  = double(nml)  / double(nc2)
   dbar_c2[i] = double(nbar) / double(nc2)
   
;;ENDFOR

myps,'./emaf_paper/plots/upsilon_vs_grs.eps',xsize=13
multiplot,[3,1]

plot,grs_dml,upsilon,ytit=cgSymbol_extra('Ypsilon'),thick=3,$
     xtit='GRS Matching Fraction for d!dML!n',xr=[0,0.99999],/xst

cgOplot,dml_c1, upsilon,thick=3,color='Dodger Blue'
cgOplot,dml_c15,upsilon,thick=3,color='Firebrick'
cgOplot,dml_c2, upsilon,thick=3,color='Forest Green'

al_legend,/top,/right,charsize=0.8,linestyle=0,box=0,thick=3,$
          color=['Opposite','Dodger Blue','Firebrick','Forest Green'],$
          ['All','C '+geq+' 0.1','C '+geq+' 0.15','C '+geq+' 0.2']

multiplot

plot,grs_dbar,upsilon,thick=3,$ ;ytit=cgSymbol_extra('Ypsilon'),$
     xtit='GRS Matching Fraction for d!dbar!n',xr=[0,0.99999],/xst

cgOplot,dbar_c1, upsilon,thick=3,color='Dodger Blue'
cgOplot,dbar_c15,upsilon,thick=3,color='Firebrick'
cgOplot,dbar_c2, upsilon,thick=3,color='Forest Green'

al_legend,/top,/right,charsize=0.8,linestyle=0,box=0,thick=3,$
          color=['Opposite','Dodger Blue','Firebrick','Forest Green'],$
          ['All','C '+geq+' 0.1','C '+geq+' 0.15','C '+geq+' 0.2']


multiplot

;; Figure out the upsilon --> R_\kappa inversion...

cgplot,charsize=1.0,[0,2500],[0,2.5],/nodata,xtitle='R!d$\kappa$!n'

vline,1000,color=cgcolor('Orange'),thick=20

colors = ['Dodger Blue','Opposite','Firebrick']
j=0
FOR Td = 15.d,25.d, 5. DO BEGIN
   Rk = upsilon * (planck_mm(299.792458d / 271.1d, Td) * 2.902d-8)
   cgoplot,Rk,upsilon,thick=3,color=colors[j]
   j++
ENDFOR

al_legend,/top,/left,color=colors,linestyle=0,box=0,thick=3,$
          ['15 K','20 K','25 K']

myps,/done,/mp

ENDFOR



END
